Simulate snow accumulation and melting
!! Simulate snow accumulation and melting !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.3 - 15th January 2025 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 15/Feb/2023 | Original code | ! | 1.1 | 29/May/2023 | change in configuration file to read initial snow water equivalent and water in snow | ! | 1.2 | 02/Dec/2024 | added refreezing of water retained in snow and snow water holding capacity | ! | 1.3 | 15/Jan/2025 | refreezing coefficient and snow water holding capacity read from configuration file | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! Module to simulate snow accumulation and melting. ! The code is a rewriting of the original code ! initially developed in the master thesis by ! Alessio Salandin and Giorgio Valè (2003) ! MODULE Snow ! Modules used: USE DataTypeSizes, ONLY : & ! Imported Type Definitions: short, float USE GridLib, ONLY : & !Imported definitions: grid_real, grid_integer, & !Imported routines: NewGrid, GridDestroy, & !Imported parameters: NET_CDF USE IniLib, ONLY : & !Imported definitions: IniList, & !Imported routines: IniOpen, IniClose, & SectionIsPresent, KeyIsPresent, & IniReadReal, IniReadInt USE LogLib, ONLY: & !Imported routines: Catch USE GridOperations, ONLY : & !Imported routines: GridByIni, CellArea, & !Imported operators: ASSIGNMENT (=) USE Chronos, ONLY : & !Imported definitions: DateTime, & !Imported operations: OPERATOR (+), & OPERATOR (==) USE AirTemperature, ONLY: & !imported variables: temperatureAir USE Precipitation, ONLY: & !imported variables: precipitationRate USE StringManipulation, ONLY : & !Imported routines: ToString USE Morphology, ONLY: & !Imported routines: DownstreamCell, CheckOutlet USE MorphologicalProperties, ONLY: & !imported variables: flowDirection, dem USE ObservationalNetworks, ONLY : & !Imported routines: ReadMetadata, CopyNetwork, & WriteMetadata, DestroyNetwork, & AssignDataFromGrid, WriteData, & !Imported defined variable: ObservationalNetwork USE Utilities, ONLY : & !Imported routines: GetUnit !Global variables: INTEGER (KIND = short) :: dtSnow !!computation time step TYPE (grid_real) :: swe !! snow water equivalent (m) TYPE (grid_real) :: rainfall !! liquid precipitation (m/s) TYPE (grid_real) :: meltCoefficient !! melt coefficient (mm/day/°C) TYPE (grid_real) :: meltTemperature !! threshold temperature for melt starts (°C) TYPE (grid_real) :: upperTemperature !! upper temperature for partitioning precipitation (°C) TYPE (grid_real) :: lowerTemperature !! lower temperature for partitioning precipitation (°C) TYPE (grid_real) :: snowConductivity !! snow hydraulic conductivity (m/s) TYPE (grid_real) :: waterInSnow !! free water in snow pack (m) TYPE (grid_real) :: swhc !!snow water holding capacity (-) TYPE (grid_real) :: cRefreeze !!coefficient of refreeze (-) TYPE (grid_real) :: excessWaterSnowFlux !! excess of water retained in snow flux (m/s) TYPE (grid_real) :: snowMelt !! snow melt amount within the time step (m) TYPE (grid_real) :: rainfallrate !! precipitation liquid fraction (m/s) !Global routines: PUBLIC :: SnowInit PUBLIC :: SnowUpdate PUBLIC :: SnowPointInit PUBLIC :: DegreeDay !local variables: TYPE (grid_integer) , PRIVATE :: meltModel !!melt model map TYPE (grid_real) , PRIVATE :: QinSnow, QoutSnow TYPE (ObservationalNetwork), PRIVATE :: sites TYPE (ObservationalNetwork), PRIVATE :: sitesSWE INTEGER (KIND = short) , PRIVATE :: fileUnitPointSWE TYPE (DateTime) , PRIVATE :: timePointExport !local routines PRIVATE :: SnowParameterUpdate PRIVATE :: SnowPointExport PRIVATE :: Refreezing !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! Initialize snow model SUBROUTINE SnowInit & ! (inifile, mask, time) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: inifile !!stores configuration information TYPE (grid_integer), INTENT(IN) :: mask TYPE (DateTime), INTENT(IN) :: time !local declarations: TYPE (IniList) :: iniDB REAL (KIND = float) :: scalar INTEGER (KIND = short) :: iscalar !---------------------end of declarations-------------------------------------- !open and read configuration file CALL IniOpen (inifile, iniDB) !load melt model IF (SectionIsPresent('melt-model', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'melt-model') ) THEN iscalar = IniReadInt ('scalar', iniDB, 'melt-model') CALL NewGrid (meltModel, mask, iscalar) ELSE CALL GridByIni (iniDB, meltModel, section = 'melt-model') END IF ELSE CALL Catch ('error', 'Snow', & 'melt-model not found in configuration file' ) END IF !set melt coefficient IF (SectionIsPresent('melt-coefficient', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'melt-coefficient') ) THEN scalar = IniReadReal ('scalar', iniDB, 'melt-coefficient') CALL NewGrid ( meltCoefficient, mask, scalar ) ELSE CALL GridByIni (iniDB, meltCoefficient, & section = 'melt-coefficient', & time = time ) END IF ELSE CALL Catch ('error', 'Snow', & 'melt-coefficient not found in configuration file' ) END IF !set threshold temperature for melt starts (°C) IF (SectionIsPresent('melt-threshold-temperature', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'melt-threshold-temperature') ) THEN scalar = IniReadReal ('scalar', iniDB, 'melt-threshold-temperature') CALL NewGrid ( meltTemperature, mask, scalar ) ELSE CALL GridByIni (iniDB, meltTemperature, & section = 'melt-threshold-temperature', & time = time ) END IF ELSE CALL Catch ('error', 'Snow', & 'melt-threshold-temperature not found in configuration file' ) END IF !set upper temperature for partitioning precipitation (°C) IF (SectionIsPresent('partitioning-upper-temperature', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'partitioning-upper-temperature') ) THEN scalar = IniReadReal ('scalar', iniDB, 'partitioning-upper-temperature') CALL NewGrid ( upperTemperature, mask, scalar ) ELSE CALL GridByIni (iniDB, upperTemperature, & section = 'partitioning-upper-temperature', & time = time ) END IF ELSE CALL Catch ('error', 'Snow', & 'partitioning-upper-temperature not found in configuration file' ) END IF !set lower temperature for partitioning precipitation (°C) IF (SectionIsPresent('partitioning-lower-temperature', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'partitioning-lower-temperature') ) THEN scalar = IniReadReal ('scalar', iniDB, 'partitioning-lower-temperature') CALL NewGrid( lowerTemperature, mask, scalar ) ELSE CALL GridByIni (iniDB, lowerTemperature, & section = 'partitioning-lower-temperature', & time = time ) END IF ELSE CALL Catch ('error', 'Snow', & 'partitioning-lower-temperature not found in configuration file' ) END IF !set snow hydraulic conductivity (m/s) IF (SectionIsPresent('hydraulic-conductivity', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'hydraulic-conductivity') ) THEN scalar = IniReadReal ('scalar', iniDB, 'hydraulic-conductivity') CALL NewGrid ( snowConductivity, mask, scalar ) ELSE CALL GridByIni ( iniDB, snowConductivity, & section = 'hydraulic-conductivity', & time = time ) END IF ELSE CALL Catch ('error', 'Snow', & 'hydraulic-conductivity not found in configuration file' ) END IF !set refreezing coefficient (-) IF (SectionIsPresent('refreezing-coefficient', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'refreezing-coefficient') ) THEN scalar = IniReadReal ('scalar', iniDB, 'refreezing-coefficient') CALL NewGrid ( cRefreeze, mask, scalar ) ELSE CALL GridByIni ( iniDB, cRefreeze, & section = 'refreezing-coefficient' ) END IF ELSE !set default value CALL NewGrid ( cRefreeze, mask, 0.05 ) END IF !set snow water holding capacity (-) IF (SectionIsPresent('water-holding-capacity', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'water-holding-capacity') ) THEN scalar = IniReadReal ('scalar', iniDB, 'water-holding-capacity') CALL NewGrid ( swhc, mask, scalar ) ELSE CALL GridByIni ( iniDB, swhc, & section = 'water-holding-capacity' ) END IF ELSE !set default value CALL NewGrid ( swhc, mask, 0.1 ) END IF !set initial optional variables !swe IF (SectionIsPresent('snow-water-equivalent', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'snow-water-equivalent') ) THEN scalar = IniReadReal ('scalar', iniDB, 'snow-water-equivalent') CALL NewGrid ( swe, mask, scalar ) ELSE CALL GridByIni (iniDB, swe, section = 'snow-water-equivalent') END IF ELSE !set to default = 0 CALL NewGrid ( swe, mask, 0. ) END IF !water in snow IF (SectionIsPresent('water-in-snow', iniDB)) THEN IF (KeyIsPresent ('scalar', iniDB, 'water-in-snow') ) THEN scalar = IniReadReal ('scalar', iniDB, 'water-in-snow') CALL NewGrid ( waterInSnow, mask, scalar ) ELSE CALL GridByIni (iniDB, waterInSnow, section = 'water-in-snow') END IF ELSE !set to default = 0 CALL NewGrid ( waterInSnow, mask, 0. ) END IF !Configuration terminated. Deallocate ini database CALL IniClose (iniDB) !allocate variables CALL NewGrid ( excessWaterSnowFlux, mask, 0. ) CALL NewGrid ( snowMelt, mask, 0. ) CALL NewGrid ( QinSnow, mask, 0. ) CALL NewGrid ( QoutSnow, mask, 0. ) RETURN END SUBROUTINE SnowInit !============================================================================== !| Description: ! compute accumulation and melting and update snow water equivalent SUBROUTINE SnowUpdate & ! ( time, mask ) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT (IN) :: time TYPE (grid_integer), INTENT (IN) :: mask !local declarations: INTEGER (KIND = short) :: i, j, is, js REAL (KIND = float) :: meltRate !!snow melt rate (m/s) REAL (KIND = float) :: refreezingRate !!refreezing rate of water retained in snowpack (m/s) !REAL (KIND = float) :: cRefreeze = 0.05!!coefficient of refreeze (-) !REAL (KIND = float) :: swhc = 0.1 !!snow water holding capacity (-) REAL (KIND = float) :: alpha !!precipitation partitioning factor REAL (KIND = float) :: liquid !!liquid fraction of precipitation (rain) (m/s) REAL (KIND = float) :: solid !!solid fraction of precipitation (snow) (m/s) REAL (KIND = float) :: tlower REAL (KIND = float) :: tupper REAL (KIND = float) :: tAir REAL (KIND = float) :: cmelt !!snow melt coefficient (m/s/°C) REAL (KIND = float) :: tmelt REAL (KIND = float) :: swe_tminus1 !!swe at previous time step REAL (KIND = float) :: ddx !!actual cell length (m) REAL (KIND = float) :: ds !!thickness of the saturated depth (m) REAL (KIND = float) :: cellWidth !!cell width (m) REAL (KIND = float) :: ijSlope !!local slope (m/m) REAL (KIND = float) :: qin, qout !!input and output water discharge in snowpack REAL (KIND = float) :: area !!cell area (m2) !------------------------------end of declarations----------------------------- !check if parameter update is required CALL SnowParameterUpdate (time) !initialize melt grid snowMelt = 0. !computer inter-cell lateral water flux QinSnow = 0. QoutSnow = 0. DO i = 1, mask % idim DO j = 1, mask % jdim IF ( mask % mat (i,j) == mask % nodata ) CYCLE !set downstream cell (is,js) CALL DownstreamCell ( i, j, flowDirection % mat(i,j), is, js, ddx, & flowDirection ) IF ( CheckOutlet (i, j, is, js, flowDirection) ) CYCLE !thickness of the saturated depth ds = waterInSnow % mat (i,j) !cell width cellWidth = ( CellArea (mask, i, j) ) ** 0.5 !local slope ijSlope = ( dem % mat (i,j) - dem % mat (is,js) ) / ddx IF ( ijSlope <= 0. ) THEN ijSlope = 0.0001 END IF QoutSnow % mat (i,j) = cellWidth * ds * & snowConductivity % mat (i,j) * ijSlope !output becomes input of downstream cell QinSnow % mat (is,js) = QinSnow % mat (is,js) + QoutSnow % mat (i,j) END DO END DO !loop across cells DO i = 1, mask % idim DO j = 1, mask % jdim !skip nodata IF ( mask % mat (i,j) == mask % nodata ) THEN CYCLE END IF !potential snow melt rate tAir = temperatureAir % mat (i,j) cmelt = meltCoefficient % mat (i,j) / 1000. / 86400. tmelt = meltTemperature % mat (i,j) SELECT CASE ( meltModel % mat (i,j) ) CASE (1) !degree-day temperature based meltRate = DegreeDay ( tAir, tmelt, cmelt ) CASE DEFAULT CALL Catch ('error', 'Snow', 'melt model not implemented: ', & argument = ToString (meltModel % mat (i,j) ) ) END SELECT ! actual snow melt, limited by available snow ! in the previous time step meltRate = MIN ( meltRate, swe % mat (i,j) / dtSnow ) snowMelt % mat (i,j) = meltRate * dtSnow !precipitation partitioning tupper = upperTemperature % mat (i,j) tlower = lowerTemperature % mat (i,j) IF ( tAir <= tlower ) THEN alpha = 0. ELSE IF ( tAir > tlower .AND. tAir < tupper ) THEN alpha = ( tAir - tlower ) / ( tupper - tlower ) ELSE alpha = 1. END IF liquid = alpha * precipitationRate % mat (i,j) solid = ( 1. - alpha ) * precipitationRate % mat (i,j) ! potential refreezing rate refreezingRate = Refreezing ( tAir, tmelt, cmelt, cRefreeze % mat (i,j) ) ! actual refreezing, limited by available retained water ! in the previous time step refreezingRate = MIN ( refreezingRate, & waterInSnow % mat (i,j) / dtSnow ) !update snow water equivalent swe % mat (i,j) = swe % mat (i,j) + & ( solid - meltRate + refreezingRate ) * dtSnow !update water in snow area = CellArea (mask, i, j) qin = QinSnow % mat (i,j) / area qout = QoutSnow % mat (i,j) / area IF ( swe % mat (i,j) > 0. ) THEN waterInSnow % mat (i,j) = waterInSnow % mat (i,j) + & snowMelt % mat (i,j) + & liquid * dtSnow + & ( qin - qout ) * dtSnow + & - refreezingRate * dtSnow rainfallRate % mat (i,j) = 0. ELSE !rainfall does not contribute waterInSnow % mat (i,j) = waterInSnow % mat (i,j) + & snowMelt % mat (i,j) + & ( qin - qout ) * dtSnow + & - refreezingRate * dtSnow rainfallRate % mat (i,j) = liquid END IF IF ( waterInSnow % mat (i,j) < 0. ) THEN waterInSnow % mat (i,j) = 0. END IF ! check retained water does not exceed water holding capacity of snow IF ( waterInSnow % mat (i,j) > swhc % mat (i,j) * swe % mat (i,j) ) THEN !exceeding water becomes runoff excessWaterSnowFlux % mat (i,j) = ( waterInSnow % mat (i,j) - & swhc % mat (i,j) * swe % mat (i,j) ) / dtSnow !set maximum waterInSnow % mat (i,j) = swhc % mat (i,j) * swe % mat (i,j) END IF END DO END DO !write point SWE IF ( time == timePointExport ) THEN CALL SnowPointExport ( time ) END IF RETURN END SUBROUTINE SnowUpdate !============================================================================== !| Description: ! update parameter map that change in time SUBROUTINE SnowParameterUpdate & ! (time) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT (IN) :: time !local declarations: CHARACTER (LEN = 300) :: filename CHARACTER (LEN = 300) :: varname !------------------------------end of declarations----------------------------- !melt coefficient IF ( time == meltCoefficient % next_time ) THEN !destroy current grid filename = meltCoefficient % file_name varname = meltCoefficient % var_name CALL GridDestroy (meltCoefficient) !read new grid in netcdf file CALL NewGrid (meltCoefficient, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !melt temperature IF ( time == meltTemperature % next_time ) THEN !destroy current grid filename = meltTemperature % file_name varname = meltTemperature % var_name CALL GridDestroy (meltTemperature) !read new grid in netcdf file CALL NewGrid (meltTemperature, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !upper temperature IF ( time == upperTemperature % next_time ) THEN !destroy current grid filename = upperTemperature % file_name varname = upperTemperature % var_name CALL GridDestroy (upperTemperature) !read new grid in netcdf file CALL NewGrid (upperTemperature, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !lower temperature IF ( time == lowerTemperature % next_time ) THEN !destroy current grid filename = lowerTemperature % file_name varname = lowerTemperature % var_name CALL GridDestroy (lowerTemperature) !read new grid in netcdf file CALL NewGrid (lowerTemperature, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF !hydraulic conductivity IF ( time == snowConductivity % next_time ) THEN !destroy current grid filename = snowConductivity % file_name varname = snowConductivity % var_name CALL GridDestroy (snowConductivity) !read new grid in netcdf file CALL NewGrid (snowConductivity, TRIM(filename), NET_CDF, & variable = TRIM(varname), time = time) END IF RETURN END SUBROUTINE SnowParameterUpdate !============================================================================== !| Description: ! compute melt rate (m/s) using degreeday temperature based model FUNCTION DegreeDay & ! ( tempAir, tempMelt, cMelt ) & ! RESULT (melt) IMPLICIT NONE !Arguments with intent(in): REAL (KIND = float), INTENT(IN) :: tempAir !! air temperature (°C) REAL (KIND = float), INTENT(IN) :: tempMelt !! melting temperature (°C) REAL (KIND = float), INTENT(IN) :: cMelt !! melting coefficient (m/s/°C) !local declarations: REAL (KIND = float) :: melt !! melt rate (m/s) !---------------------end of declarations-------------------------------------- IF ( tempAir <= tempMelt ) THEN melt = 0. ELSE melt = cMelt * ( tempAir - tempMelt ) END IF RETURN END FUNCTION DegreeDay !============================================================================== !| Description: ! compute refreezing rate (m/s) of the water retained in snowpack ! ! References: ! van Verseveld, W. J., Weerts, A. H., Visser, M., Buitink, J., ! Imhoff, R. O., Boisgontier, H., Bouaziz, L., Eilander, D., Hegnauer, M., ! ten Velden, C., and Russell, B.: Wflow_sbm v0.7.3, a spatially distributed ! hydrological model: from global data to local applications, Geosci. ! Model Dev., 17, 3199–3234, https://doi.org/10.5194/gmd-17-3199-2024, 2024 ! FUNCTION Refreezing & ! ( tempAir, tempMelt, cMelt, cRefreeze ) & ! RESULT (freeze) IMPLICIT NONE !Arguments with intent(in): REAL (KIND = float), INTENT(IN) :: tempAir !! air temperature (°C) REAL (KIND = float), INTENT(IN) :: tempMelt !! melting temperature (°C) REAL (KIND = float), INTENT(IN) :: cMelt !! melting coefficient (m/s/°C) REAL (KIND = float), INTENT(IN) :: cRefreeze !! refreezing coefficient (-) !local declarations: REAL (KIND = float) :: freeze !! refreezing rate (m/s) !---------------------end of declarations-------------------------------------- IF ( tempAir <= tempMelt ) THEN freeze = cRefreeze * cMelt * ( tempMelt - tempAir ) ELSE freeze = 0. END IF RETURN END FUNCTION Refreezing !============================================================================== !| Description: ! Initialize export of point site data SUBROUTINE SnowPointInit & ! ( pointfile, path_out, time ) IMPLICIT NONE !Arguments with intent (in): CHARACTER (LEN = *), INTENT(IN) :: pointfile !!file containing coordinate of points CHARACTER (LEN = *), INTENT(IN) :: path_out !!path of output folder TYPE (DateTime), INTENT(IN) :: time !!start time !local declarations INTEGER (KIND = short) :: err_io INTEGER (KIND = short) :: fileunit !-------------------------end of declarations---------------------------------- timePointExport = time !open point file fileunit = GetUnit () OPEN ( unit = fileunit, file = pointfile(1:LEN_TRIM(pointfile)), & status='OLD', iostat = err_io ) IF ( err_io /= 0 ) THEN !file does not exist CALL Catch ('error', 'Snow', 'out point file does not exist') END IF !Read metadata CALL ReadMetadata (fileunit, sites) !check dt IF (.NOT. ( MOD ( sites % timeIncrement, dtSnow ) == 0 ) ) THEN CALL Catch ('error', 'Snow', 'dt out sites must be multiple of dtSnow ') END IF CLOSE ( fileunit ) !create virtual network and initialize file for output fileUnitPointSWE = GetUnit () OPEN ( unit = fileUnitPointSWE, & file = TRIM(path_out) // 'point_swe.fts' ) CALL CopyNetwork ( sites, sitesSWE ) sitesSWE % description = 'snow water equivalent data exported from FEST' sitesSWE % unit = 'mm' sitesSWE % offsetZ = 0. CALL WriteMetadata ( network = sitesSWE, & fileunit = fileUnitPointSWE ) CALL WriteData (sitesSWE, fileUnitPointSWE, .TRUE.) ! destroy sites CALL DestroyNetwork ( sites ) RETURN END SUBROUTINE SnowPointInit !============================================================================== !| Description: ! Export of point site data SUBROUTINE SnowPointExport & ! ( time ) IMPLICIT NONE !Arguments with intent(in): TYPE (DateTime), INTENT (IN) :: time !local declarations: INTEGER (KIND = short) :: i !-------------------------end of declarations---------------------------------- !set current time sitesSWE % time = time !populate data CALL AssignDataFromGrid ( swe, sitesSWE, scaleFactor = 1000. ) !write data CALL WriteData ( sitesSWE, fileUnitPointSWE ) timePointExport = timePointExport + sitesSWE % timeIncrement RETURN END SUBROUTINE SnowPointExport END MODULE Snow